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Abstract 

We  develop  in  this  article  an  improved  version  of  the  fifth-order  weighted  essentially  non-oscillatory 
(WENO)  scheme.  Through  the  novel  use  of  higher  order  information  already  present  in  the  framework 
of  the  classical  scheme,  new  smoothness  indicators  are  devised  and  we  obtain  a  new  WENO  scheme 
with  less  dissipation  than  the  classical  WENO  of  Jiang  and  Shu  [2],  with  the  same  computational  cost, 
and  a  slightly  better  performance  than  the  improved  mapped  version  of  Henrick  et  al  [3].  We  show 
that  the  enhancements  of  the  new  scheme  come  from  its  ability  to  assign  substantially  larger  weights 
to  discontinuous  stencils  than  the  previous  versions  of  WENO.  Nnmerical  experiments  with  the  linear 
advection  of  discontinuons  fnnctions  and  the  one  dimensional  Euler  system  of  equations  are  conducted  to 
demonstrate  the  benefit  of  using  this  improved  version  of  the  WENO  scheme  for  hyperbolic  conservation 
laws. 
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1  Introduction 

In  the  numerical  simulation  of  compressible  flows  modeled  by  means  of  hyperbolic  conservation  laws  in  the 

form 


—  +  V-F{u)  =  0,  (1) 

the  development  of  finite  time  discontinuities  generates  0(1)  oscillations,  known  as  the  Gibbs  phenomenon, 
causing  loss  of  accuracy  and  numerical  instability.  Among  many  choices  of  shock  capturing  numerical 
schemes  such  as  the  Piecewise  Parabolic  method  (PPM)  [7],  the  Essentially  Non-Oscillatory  scheme  (ENO) 
[6],  high  order  Weighted  Essentially  Non-Oscillatory  schemes  (WENO)  [1,  2]  have  been  extensively  used  for 
the  simulation  of  the  fine  scale  and  delicate  structures  of  the  physical  phenomena  related  to  shock-turbulence 
interactions. 

WENO  schemes  owe  their  success  to  the  use  of  a  dynamic  set  of  stencils,  where  a  nonlinear  convex 
combination  of  lower  order  polynomials  adapts  either  to  a  higher  order  approximation  at  smooth  parts 
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of  the  solution,  or  to  an  upwind  spatial  discretization  that  avoids  interpolation  across  discontinuities  and 
provides  the  necessary  dissipation  for  shock  capturing.  The  nonlinear  coefficients  of  the  convex  combination, 
hereafter  referred  to  as  weights,  are  based  on  the  local  smoothness  indicators,  which  measure  the  sum  of  the 
normalized  squares  of  the  scaled  norms  of  all  derivatives  of  the  lower  order  polynomials  [2] .  An  essentially 
zero  weight  is  assigned  to  those  lower  order  polynomials  whose  underlining  stencils  contain  high  gradients 
and/or  shocks,  yielding  an  essentially  non-oscillatory  solution  at  discontinuities.  At  smooth  regions,  higher 
order  is  achieved  through  the  mimicking  of  the  central  upwinding  scheme  of  maximum  order,  when  all 
smoothness  indicators  are  about  the  same  size.  Hence,  an  efficient  design  of  these  smoothness  indicators  is 
a  very  important  issue  for  WENO  schemes. 

The  classical  choice  of  smoothness  indicators  in  [2]  generated  weights  that  failed  to  recover  the  maximum 
order  of  the  scheme  at  critical  points  of  the  solution.  This  fact  was  clearly  pointed  out  at  Henrick  et  al. 
[3].  In  their  study,  necessary  and  sufficient  conditions  on  the  weights,  for  optimality  of  the  order,  were 
derived  and  a  correcting  mapping  to  be  applied  to  the  classical  weights  was  devised.  The  resulting  mapped 
WENO  scheme  of  [3]  recovered  the  optimal  order  of  convergence  at  critical  points  and  presented  sharper 
results  close  to  discontinuities.  In  this  article,  we  follow  a  different  approach,  which  is  to  improve  on  the 
classical  smoothness  indicators  to  obtain  weights  that  satisfies  the  sufficient  conditions  for  optimal  order. 
Taylor  series  analysis  of  the  classical  smoothness  indicators  reveals  that  a  simple  combination  of  them  would 
give  higher  order  information  about  the  regularity  of  the  numerical  solution.  The  incorporation  of  this 
new  higher  order  information  into  the  weights  definition  improves  the  convergence  order  at  smooth  parts 
of  the  solution,  as  well  as  decreases  dissipation  close  to  discontinuities,  while  maintaining  stability  and  an 
essentially  non-oscillatory  behavior. 

The  enhancements  of  the  new  scheme  come  from  the  bigger  weights  that  it  assigns  to  discontinuous 
stencils.  Contrary  to  common  belief,  the  strategy  should  be  to  augment  the  influence  of  the  stencil  containing 
the  discontinuity  as  much  as  possible,  without  destroying  the  essentially  non-oscillatory  behavior.  A  com¬ 
parison  of  the  weights  of  the  classical,  the  mapped  and  the  new  WENO  scheme  close  to  a  discontinuity  shows 
that  the  ratio  between  the  weight  of  a  discontinuous  stencil  and  a  continuous  one  increases  slightly  from  the 
classical  weights  to  the  mapped  weights,  and  increases  substantially  with  the  new  weights  proposed  in  this 
article.  This  was  made  possible  through  the  use  of  higher  order  smoothness  indicators  already  available  at 
the  definition  of  the  classical  weights. 

This  paper  is  organized  as  follows:  In  Section  2,  the  classical  WENO  scheme  of  Jiang  and  Shu  [2]  and 
the  mapped  weights  version  of  Henrick  et  al.  [3]  are  described  and  some  of  the  relevant  analytical  results 
are  reviewed.  The  new  WENO  scheme  is  introduced  at  Section  3  where  a  detailed  discussion  about  the  new 
smoothness  indicators  is  given.  In  Section  4,  we  perform  a  numerical  comparison  of  all  the  schemes  in  the 
linear  advection  of  discontinuous  functions  and  in  the  solution  of  the  one  dimensional  Euler  equations  in  the 
classical  cases  of  the  shock  density  wave  interaction  and  the  interactive  blastwaves  problems. 


2  Weighted  essentially  non-oscillatory  schemes 


In  this  section  we  describe  the  fifth  order  weighted  essentially  non-oscillatory  conservative  finite  difference 
scheme  when  applied  to  hyperbolic  conservation  laws  as  in  (1).  Without  loss  of  generality,  we  will  restrict 
our  discussion  in  this  article  to  the  one  dimensional  scalar  case.  Extensions  to  system  of  equations  and 
higher  spatial  dimensions  present  no  extra  complexity  with  regards  to  our  main  point  which  is  the  design  of 
new  smoothness  indicators  for  the  WENO  scheme. 

Consider  an  uniform  grid  defined  by  the  points  Xi  =  iAx,  i  =  0, . . . ,  N ,  which  are  also  called  cell  centers, 
with  cell  boundaries  given  by  Xj+i  =  Xi  +  The  semi-discretized  form  of  (1)  by  the  method  of  lines, 
yields  a  system  of  ordinary  differential  equations 


dui{t) 

dt 


(2) 


where  Ui{t)  is  a  numerical  approximation  to  the  point  value  u{xi,t). 
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Figure  1:  The  computational  uniform  grid  Xi  and  the  total  5-points  stencil  5®,  composed  of  the  three  3-points 
stencils  Sk,k  =  0, 1,  2,  used  for  the  fifth  order  WENO  reconstruction. 


A  conservative  finite-difference  formulation  for  hyperbolic  conservation  laws  requires  high-order  consis¬ 
tent  numerical  fluxes  at  the  cell  boundaries  in  order  to  form  the  flux  difference  across  the  uniformly-spaced 
cells.  The  conservative  property  of  the  spatial  discretization  is  obtained  by  implicitly  defining  the  numerical 
flux  function  h{x)  as 


1 

Ax 


such  that  the  spatial  derivative  in  (2)  is  exactly  approximated  by  a  conservative  finite  difference  formula  at 
the  cell  boundaries, 


dui{t) 

dt 


—  ( 
Ax  V 


-  K 


(3) 


where  =  h{x  ±  ^). 

High  order  polynomial  interpolations  to  are  computed  using  known  grid  values  of  /.  The  classical 
fifth-order  WENO  scheme  uses  a  5-points  stencil,  hereafter  named  S'®,  which  is  subdivided  into  three  other 
3-points  stencils  So,Si,S2,  as  shown  in  Fig.  1.  The  fifth-order  polynomial  approximation  "Ci+i  =  ^i+i  + 

0(Ax®)  is  built  through  the  convex  combination  of  interpolated  values  fk  at  by  the  third- 

degree  polynomials  /fc,  defined  in  each  one  of  the  stencils  Sk- 


2 

Vi+I  =  '^UJkfk  > 


(4) 


where 


2 

fk  Ckjfi-k+j 

j=o 


=  h  +0{Ax^). 


(5) 


The  Ckj  are  Lagrangian  interpolation  coefficients  (see  [2]),  which  depend  on  the  left-shift  parameter  fc  =  0, 1,  2, 
but  not  on  the  values  fi.  The  weights  ujk  are  defined  as 


Oik 


dk 

(/3fe  +  e)P  ■ 


(6) 
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We  refer  to  as  the  unnormalized  weights.  The  parameter  e  is  used  to  avoid  the  division  by  zero  in  the 
denominator  and  p  =  2  is  chosen  to  increase  the  difference  of  scales  of  distinct  weights  at  non-smooth  parts 
of  the  solution. 

The  smoothness  indicators  f3k  measure  the  regularity  of  the  fc_th  polynomial  approximation  fk{x)  at 


the  Stencil  Sk  and  are  given  by 

/ -1  ^  X  .  1  \  / 

l  —  ±  i  ^ 

The  expression  of  the  fUk  in  terms  of  the  cell  values  of  /  are  given  by 

Po  =  g(/.-2-2/,_i+/,)V  J(/i_2-4/,_i+3/.f ,  (8) 

Pi  =  -^{fi-l  —  +  ft+l)"^  +  ^ifi-l  —  fi+l)"^  :  (9) 

P2  =  +  fi+2)^  +  ^{3fi  —  4:fi+l  +  fi+2f  ■  (10) 


The  coefficients  do  =  di  =  |,  ^2  =  in  (6)  are  called  the  ideal  weights  for  the  convex  combination 
(4)  since  they  generate  the  central  upstream  fifth-order  scheme  for  the  5-points  stencil  S^.  The  general  idea 
of  the  weights  definition  (6)  is  that  on  smooth  parts  of  the  solution  the  smoothness  indicators  Pk  are  all 
small  and  about  the  same  size,  generating  weigths  ojk  that  are  good  approximations  to  the  ideal  weights 
dfc.  On  the  other  hand,  if  the  stencil  Su  contains  a  discontinuity,  Pk  is  0(1)  and  the  corresponding  weight 
LOk  is  small  relatively  to  the  other  weights.  This  implies  that  the  influence  of  the  polynomial  approximation 
of  taken  across  the  discontinuity  is  diminished  up  to  the  point  where  the  convex  combination  (4)  is 
essentially  non-oscillatory.  Figure  1  shows  the  case  where  stencil  S'2  is  discontinuous,  yielding  /3o  and  Pi  to 
be  much  smaller  than  P2.  By  (6),  this  results  on  ui^  being  a  small  number  in  the  convex  combination  (4) 
(See  also  Fig.  3(a)  of  Section  3). 

The  process  described  above  is  called  the  WENO  reconstruction  step,  for  it  reconstructs  the  values  of 
h{x)  at  the  cell  boundaries  of  the  interval  h  =  from  its  cell  averaged  values  f{x)  in  the  intervals 

S'fe,  /c  =  0, 1,  2.  In  [3],  truncation  error  analysis  of  the  finite  difference  equation  (3)  led  to  sufficient  conditions 
on  the  weights  uik  for  the  WENO  scheme  to  achieve  the  formal  fifth-order  of  convergence  at  smooth  parts  of 
the  solution.  It  was  found  that  at  critical  points,  points  where  the  hrst  derivative  of  the  solution  vanishes, 
convergence  degraded  to  only  third  order,  a  fact  that  was  being  hidden  by  the  uniformization  of  the  weights 
caused  by  the  use  of  a  large  value  for  e  in  (6).  Since  this  analysis  is  relevant  to  the  description  of  the  new 
WENO  scheme  introduced  at  next  section,  we  recall  below  its  most  important  steps. 

Substituting  the  convex  combination  (4)  into  the  finite  difference  formula  (3)  we  obtain: 


v,+  i  -v,_i 


k^O  k^O 

/  2  2  \  2  2 

=  dkV,_i  +  -  dk)vi+i  -  Y^^k  -  dk)vi_i 


=  -dk) 
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~Y^^^k  -  dk)  hj_i-l-BfcAx" 


=  v'{xi)/S.x  -\-  0{Ax^)  +  y^(i^fc  -  dk)AkAx^  —  YX^^k  -dk)BkAx^. 


fe=0 


fc=0 


Hence,  the  trunction  error  becomes 


1 

Ax 


Z  Z 

+  O(Ax^)  +  -  dk)AkAx‘^  -  -  dk)BkAx‘^,  (11) 


4 


where  v^j_i  are  the  fifth-order  central  upstream  approximations  to  and  and  are  the  WENO 
weights  at  and  x^_i,  respectively,  Ak  and  are  constants  independent  of  Ax.  It  should  be  pointed 

2  2 

out  that  due  to  normalization  ^  =  '^  dk  and  that  O(Ax^)  is  still  obtained  after  the  division  by  Ax  in 

fc=0  fc=o 

the  difference  of  the  central  approximations  ,  because  of  cancellation  of  identical  0(Ax®)  terms.  Thus, 
a  sufficient  condition  on  the  weights  (6)  to  get  a  fifth-order  truncation  error  at  (3)  is  given  by 

u;±  =  dk  +  0{Ax^),  fc  =  0,l,2.  (12) 

Considering,  for  simplicity,  that  e  =  0,  it  is  easily  seen  from  (6)  that  if  the  smoothness  indicators  satisfy 
f3k  =  D{1  +  0(Ax^)),  for  a  constant  D  independent  of  fc,  then  the  unnormalized  weights  also  satisfy  au  = 
D{\  +  0(Ax^))  and  u!k  =  dk  +  O(Ax^)  (see  [3]). 

It  is  easily  seen  that  for  the  classical  scheme,  and  a  sufficiently  smooth  solution, 

Pk  =  D{1  +  0(Ax^))  and  LOk  =  dk  +  O(Ax^),  (13) 


which  does  not  satisfy  (12),  yielding  convergence  order  less  than  5.  In  [3],  the  weights  LUk  were  modified 
by  a  mapping  function  that  increased  the  approximation  to  the  ideal  weights  dk  to  the  required  third  order 
O(Ax^).  The  mapping  function  gk{i^)  was  defined  as 


,  ,  oj  (dk  +  dl  —  3dkio  + 

=  ,2  ^ 

dl+uj(l  -  2dk) 

It  is  a  nondecreasing  monotone  function  with  the  following  properties: 

•  0  <  gk(uj)  <  1,  gk(0)  =  0  and  gk(l)  =  1. 

•  gk{oj)  w  0  if  w  «  0;  gk{x>)  «  1  if  w  w  1. 

•  gk{dk)  =  dk,  g'kidk)  =  g'kidk)  =  0. 

•  gk{oj)  =  dk  +  0(Ax^’'),  if  w  =  dfc  +  0{Ax^). 


(14) 


At  critical  points  of  the  solution,  the  order  of  the  classical  scheme  decreases  more  since  we  only  obtain 
Pk  =  D{1  +  0(Ax)),  if  f"  7^  0,  and  Pk  =  D{1  +  0(1)),  if  /"  =  0  (see  Equations  (15)-(17)  below).  Thus, 
the  mapping  does  improve  the  order  of  the  classical  scheme  at  the  critical  points  of  first  order,  but  it  cannot 
help  when  too  many  derivatives  vanish. 

For  problems  with  shocks,  the  0(1)  truncation  error  at  the  discontinuities  diminishes  the  advantages  of 
such  order  improvements.  Nevertheless,  the  results  obtained  by  the  mapped  scheme  are  clearly  superior  to 
the  classical  scheme  results.  Distinctly  from  [3],  we  do  not  credit  the  enhancements  of  the  numerical  solution 
to  the  higher  order  of  the  mapped  scheme  at  critical  points,  but  to  the  smaller  dissipation  that  it  provides 
when  assigning  larger  weights  to  discontinuous  stencils,  as  we  shall  see  in  the  next  section.  We  will  also  see 
that  the  new  smoothness  indicators  proposed  in  this  work  yield  even  larger  weights  than  the  mapped  WENO 
to  stencils  with  discontinuities,  generating  even  better  solutions. 


3  An  Improved  WENO  Scheme 

In  this  section  we  devise  a  new  set  of  weights  that  satisfies  the  sufficient  condition  (12)  and  recovers  the 
formal  order  of  accuracy  at  the  smooth  regions  of  the  solution.  The  novel  idea  is  the  use  of  higher  order 
smoothness  indicators  at  the  formula  for  the  weights  of  the  convex  combination  (4).  In  order  to  define  the 
new  smoothness  indicators  of  higher  order,  we  need  the  following  Taylor  series  expansions  for  the  smoothness 
indicators  Pk  of  (7): 
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Ax®  +  0(Ax®), 


,//2  I  ^ 


■*■  0{Ax°), 
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/92  =  x'"Ax"+  Ax^+  ^x"x"'--xV'"  Ax®  +  0(Ax®). 
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The  higher  order  new  smoothness  indicator  T5  is  defined  as  the  difference  between  /Sq  and  (32,  namely 


(15) 

(16) 
(17) 


T5  =  max(|/3o  -  /32|  ,e) , 


(18) 


where  e  is  a  small  number  used  to  ensure  that  T5  ^  0.  We  see  from  (10)  and  from  equations  (15)-(17)  that 
T5  is  a  measure  of  the  higher  derivatives  of  the  function,  when  they  exist,  and  is  computed  using  the  whole 
5-points  stencil  5®.  The  following  properties  of  are  important  for  the  definition  of  the  new  WENO  scheme: 

•  If  5'®  contains  neither  discontinuities,  nor  critical  points,  then  =  0(Ax®)  <^C  (3^,  for  fc  =  0, 1,  2; 

•  if  the  solution  is  continuous  at  some  of  the  Sk,  but  discontinuous  in  the  whole  S'®,  then  /3fc  <C  T5,  for 
those  k  where  the  solution  is  continuous; 

•  T5  <  maxfe  Pk- 


Figure  2  shows  the  values  of  Pk  and  for  the  discontinuous  function: 


M(a:,0)  =  /(x) 


—  sin(7rx)  —  ^x®  — 1  <  x  <  0 

—  sin(7rx)  —  ^x®  -f  1  0  <  x  <  1 


(19) 


consisting  of  a  piecewise  sine  function  with  a  jump  discontinuity  at  x  =  0.  Note  that  ts  is  only  comparable 
to  one  of  the  Pk  at  stencils  that  include  the  discontinuity. 


X 

Figure  2:  Values  of  Pk  and  T5  for  the  discontinuous  periodic  function  (19). 

The  following  notation  will  be  used  in  order  to  distinguish  between  the  3  different  WENO  schemes 
considered  in  this  work.  The  mapped  weights  of  [3]  are  denoted  as  uj^  and  the  resulting  mapped  WENO 
scheme  as  WENO-M.  The  new  WENO  scheme,  introduced  below,  is  referred  as  WENO-Z  and  the  superscript 
z  is  added  to  all  quantities  related  to  it.  To  keep  a  coherent  notation  throughout  the  article,  the  classical 
weights  and  smoothness  indicators  carry  no  superscript,  although,  the  classical  WENO  scheme  of  [2]  is 
referred  as  WENO-JS. 
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The  new  smoothness  indicators  are  defined  as 


/3fc  +  e 

+  T5  ’ 


and  the  new  weights  as 


Wfc  = 


z^;=o 


=  4>  *  =  0,1,2. 
I^k 


(20) 


(21) 


Rigorously  speaking,  are  not  smoothness  indicators,  for  they  are  all  close  to  1  at  smooth  parts  of  the 
solution.  They  are  in  fact  the  normalization  of  the  classical  smoothness  indicators  by  higher  order  information 
contained  in  the  larger  stencil  S^. 


Remark  1  Note  thatp  =  1  in  the  definition  of  the  unnormalized  weights  a^.  Numerical  experiments  showed 
that  the  use  of  p  =  2,  as  in  the  definition  of  the  classical  unnormalized  weights  oik,  led  to  an  unstable  scheme, 
for  it  amplified  too  much  the  ratio  (24).  This  is  the  opposite  way  that  was  taken  in  the  definition  of  the 
classical  weights  where  p  =  2  was  used  to  increase  the  distance  between  a  continuous  and  a  discontinuous 
stencil  and,  therefore,  decrease  the  importance  of  this  last  one  at  the  final  convex  combination. 


Remark  2  The  role  of  e  was  discussed  in  [3],  where  it  was  shown  that  the  large  value  o/10“®,  commonly 
suggested  in  the  literature,  would  dominate  over  the  smoothness  indicators  fik  at  the  denominators  of  the 
definition  of  the  classical  weights.  Due  to  (13),  a  smaller  value  for  e  degrades  the  order  of  the  classical 
scheme,  particularly  at  critical  points.  In  this  work,  we  use  much  smaller  values  of  e  for  all  schemes,  in 
order  to  make  apparent  their  real  truncation  errors.  In  other  words,  the  parameter  it  is  reduced  to  play  only 
its  original  role  of  not  allowing  vanishing  denominators  at  the  weights  definitions. 

In  the  following,  we  show  that  Lof  =  dk  +  0{Ax^)  at  smooth  parts  of  the  solution  and  that  bigger 
weights  are  assigned  to  discontinuous  stencils  with  the  new  scheme  WENO-Z.  We  take  e  =  0  at  the  analyses 
below  due  to  the  reasons  given  in  Remark  2. 

We  now  compute  the  order  of  approximation  of  the  weights  in  (21)  to  the  ideal  weights  dk.  For  the  sake 
of  clarity,  we  will  not  consider  the  case  of  critical  points  initially.  Thus,  at  smooth  regions  of  the  solution, 
not  containing  critical  points,  we  see  from  Equations  (15)-(17)  that  fik  ~  0{Ax^),  yielding 


It  follows  that 


Wfe 


=  {l  +  T5/3fi^)  ^  =  l  +  0{Ax^),  k  =  0,1,2. 


dk  (l  +  T5/3fc 

*^1  (l  + '^5/3;  ) 


4(l  +  Q(Ax^)) 
Ei=o  'ii  (1  +  0(Aa;^)) 


dk  +  0{Ax^),  k  =  0,1,2. 


(22) 


(23) 


Thus,  the  new  weights  uf  satisfy  the  sufficient  condition  (12),  providing  the  formal  fifth  order  of  accuracy 
to  the  WENO-Z  scheme  at  the  smooth  regions  of  the  solution. 

Next,  we  examine  the  behavior  of  cof  on  stencils  containing  discontinuities  with  respect  to  the  classical 
weights.  The  analysis  can  be  performed  by  looking  at  the  behavior  of  the  smoothness  indicators  /3f.  Consider 
the  simple  case  of  a  shock  localized  in  stencil  S2,  while  the  solution  in  stencils  Sq  and  S'!  are  smooth  (see  Fig. 
1),  the  ratios  between  the  smaller  smoothness  indicators  fik,  *  =  0,1  and  the  largest  one,  (32,  are  increased 
with  for  the  new  smoothness  indicators  : 


(3f  (3k  (32  +T5  ^  (3k 

/3|  (32  (3k +  T5-  (32^ 


(24) 


using  the  fact  that  (32  >  (3k,  k  =  0,1. 
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Similarly,  the  relative  importance  of  the  stencil  S2  containing  the  discontinuity  in  the  new  convex 
combination  is  also  increased,  since: 


.,2 


dkj3? 


^  ,  fc  =  0,l, 


(25) 


Let  us  now  compare  the  sizes  of  the  weights  ojk  for  the  classical  WENO-JS,  WENO-M  and  the  new 
WENO-Z  schemes  at  the  first  step  of  the  numerical  solution  of  the  wave  equation  ut  =  Ux,  with  periodic 
initial  condition  given  by  (19). 


Eigure  3:  The  distribution  of  the  ideal  weights  dk  and  the  weights  ujk,k  =  0,1,2  for  (a)  WENO-JS,  (b) 
WENO-M  and  (c)  WENO-Z  schemes  at  the  first  step  of  the  numerical  solution  of  the  wave  equation  ut  =  Ux, 
with  periodic  initial  condition  given  by  (19).  The  ideal  weights  d/c,  fc  =  0, 1,  2  are  shown  with  line  styles  and 
the  weights  LOk,  fc  =  0, 1,  2  are  shown  with  symbols.  The  vertical  axes  are  shown  in  logj^Q  scale. 

Figures  3(a)-(c)  show  the  weights  ujk,  for  the  WENO-JS,  WENO-M  and  WENO-Z  schemes  plotted 
together  with  the  ideal  weights  d^, where  the  vertical  axis  is  in  log^g  scale.  Away  from  the  discontinuity,  at 
a;  =  0,  the  weights  ujk  (symbols)  for  all  schemes  correctly  match  the  corresponding  ideal  weights  dk  (lines). 
The  first  location  where  the  five-points  stencil  encounter  the  discontinuity  is  at  a:  =  —0.01.  At  this  grid 
point,  the  rightmost  stencil  S2  has  its  weight  0J2  decreased  to  a  much  smaller  value  than  wq  and  oji,  while 
these  two  are  slightly  increased  to  reflect  their  larger  relevance  at  the  reconstruction  step.  At  x  =  —0.005, 
the  discontinuity  is  present  at  Si  and  S2  and  a  small  value  is  also  assigned  to  oji,  yielding  wq  «  0(1)-  A 
symnmetric  scenario  occurs  at  a;  =  0.005,  where  UJ2  assumes  the  largest  value.  While  this  situation  is  general 
for  all  schemes,  the  main  difference  is  at  the  ratios  between  the  weights  for  discontinuous  and  continuous 
stencils,  as  discussed  above.  While  WENO-JS  sets  very  small  values  for  the  discontinuous  stencils,  around 
10“®,  and  the  mapping  of  WENO-M  generates  a  small  increase  on  these  values,  WENO-Z  yields  a  more 
substantial  increase  to  10“^. 

Figure  4  shows  the  numerical  solutions  of  the  wave  equation  ut  =  Ux  &t  t  =  8,  for  all  three  schemes,  along 
with  the  exact  solution.  Note  that  the  numerical  solution  generated  by  WENO-Z  is  qualitatively  similar  to 
the  one  generated  by  WENO-M,  being  slightly  sharper  at  the  discontinuity.  It  can  also  be  observed  that 
both  schemes  show  distinctly  better  results  than  the  classical  WENO-JS. 

4  Numerical  Experiments 

In  this  section  we  compare  the  numerical  performance  of  the  new  weno  scheme,  WENO-Z,  with  the  classical 
WENO-JS  and  its  improved  version  WENO-M.  We  show  at  first  that  WENO-Z  is  indeed  a  fifth-order 
scheme  by  showing  its  errors  with  a  smooth  scalar  advection  problem.  We  also  compare  its  behavior  in 


Figure  4:  Numerical  solutions  of  the  linear  wave  equation  with  the  discontinuous  initial  condition  (19)  at 
t  =  8  for  all  three  schemes.  The  exact  solution  is  shown  with  a  solid  line. 


the  advection  of  discontinuous  functions  with  the  other  two  schemes.  Finally,  the  one  dimensional  Euler 
equations  are  solved  for  the  Mach  3  shock  density  wave  interaction  and  the  interactive  blastwaves  problems. 

4.1  Linear  advection  problems 

In  this  section,  the  new  scheme,  WENO-Z,  is  used  to  solve  the  linear  transport  equation  with  smooth  and 
discontinuous  initials  conditions.  Table  I  shows  the  Li,  L2  and  Loo  errors,  along  with  the  respective  orders 
of  convergence,  when  WENO-Z  is  applied  to  the  numerical  solution  of  the  linear  advection  of  the  scalar 
smooth  function: 


u{x,  0)  =  sin(27ra:). 

One  can  see  that  WENO-Z,  with  the  new  weights  achieves  fifth-order  convergence  for  smooth  problems. 


N 

Li  Error 

Order 

L2  Error 

Order 

Loo  Error 

Order 

25 

7.9e-l 

6.2e-l 

6.5e-l 

50 

3.5e-2 

4.50 

2.8e-2 

4.49 

2.8e-2 

4.55 

100 

l.le-3 

4.98 

8.8e-4 

4.98 

8.8e-4 

4.98 

200 

3.5e-5 

5.00 

2.7e-5 

5.00 

2.7e-5 

5.00 

400 

l.le-6 

5.00 

8.6e-7 

5.00 

8.6e-7 

5.00 

Table  I:  The  Li,  L2  and  Loo  errors  of  the  linear  wave  equation,  with  a  smooth  initial  condition,  as  computed 
by  the  WENO-Z  scheme  with  the  new  weights 

Next,  WENO-Z  is  tested  with  the  linear  transport  of  discontinuous  functions,  in  the  case  of  an  initial 
condition  consisting  of  a  Gaussian,  a  triangle,  a  square-wave  and  a  semi-ellipse,  given  by 


i  [G{x,  l3,z  -  6)  +  4:G{x,  13,  z)  F  G{x,  l3,z  +  d)] 

,  x  S  [—0.8,  —0.6] 

1 

,x&  [-0.4, -0.2] 

u(x,  0)  = 

< 

1  -  |10(a:- 0.1)1 

,a;e  [0,0.2] 

(26) 

1  [7^(0;,  a,a  —  S)  +  4:F(x,  a,  a)  +  F{x,  a,  a  -I-  d)] 

,  X  S  [0.4, 0.6] 

0 

V. 

,  otherwise 

G{x,P,z)  = 

e 

5 

(27) 

F(x,a,a)  = 

•\/max(l  —  a‘^{x  —  a)2,0) , 

(28) 
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(a)  Solution 


(b)  Error 


X 


Figure  5:  Numerical  solution  and  absolute  pointwise  error  of  the  advection  equation  with  the  discontinuous 
initial  condition  (28)  as  computed  by  WENO-JS,  WENO-M  and  WENO-Z  with  N  =  200  at  t  =  8. 


(a)  Solution 


(b)  Error 


Figure  6:  Numerical  solution  and  absolute  pointwise  error  of  the  advection  equation  with  the  discontinuous 
initial  condition  (28)  as  computed  by  WENO-JS,  WENO-M  and  WENO-Z  with  N  =  400  at  t  =  8. 
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where  the  constants  are  z  =  —0.7,5  =  0.005, /3  =  ^7,0  =  0.5  and  a  =  10. 

Figures  5  and  6  show  that  WENO-Z  behaves  quantitatively  and  qualitatively  equivalent  to  WENO- 
M  with  regards  to  the  improvements  over  WENO-JS,  but  WENO-Z  still  shows  the  best  results  of  all  three 
schemes.  Also,  as  indicated  in  the  figures,  the  smaller  dissipation  of  WENO-Z  generates  slightly  better  results 
on  the  various  corners  and  pikes  of  the  solution.  Note  also  that  the  theoretical  deficiency  of  WENO-Z  at  the 
critical  points  does  not  seem  to  influence  the  performance  at  the  capturing  of  the  several  discontinuities  of 
the  solution.  This  is  expected  since  the  mapping  of  WENO-M  increases  the  order  from  3  to  5  at  the  critical 
points  for  smooth  solutions  only,  while  at  the  discontinuities,  the  overall  global  rate  of  convergence  of  the 
WENO  schemes  is  only  first  order. 


4.2  The  Euler  equations 

In  this  section  we  present  numerical  experiments  with  the  system  of  Euler  Equations  for  gas  dynamics  in 
strong  conservation  form: 

Qt  +  Fa;  =  0,  (29) 

where 

Q  =  (p,pu,E)",  F  =  (pu,pu"  +  F,(^;  +  F)u)^  (30) 

and  the  equation  of  state  is  given  by 

P=(7-l)(^E+ipu2^,  7=1.4.  (31) 

The  p,  u,  P,  E  are  the  density,  velocity,  pressure  and  total  energy  respectively.  Following  [1],  the  hyperbolicity 
of  the  Euler  equations  admits  a  complete  set  of  right  and  left  eigenvectors  for  the  Jacobian  of  the  system. 
The  eigenvalues  and  eigenvectors  are  obtained  via  the  linearized  Riemann  solver  of  Roe  [4]  and  the  first 
order  Lax-Eriedrichs  flux  is  used  as  the  low  order  building  block  for  the  high  order  reconstruction  step 
of  WENO  (see  equation  (2.5)  in  [1]).  After  projecting  the  fluxes  on  the  characteristic  fields  via  the  left 
eigenvectors,  the  high  order  WENO  reconstruction  step  is  applied  to  obtain  the  high  order  approximation 
at  the  cell  boundaries,  which  are  projected  back  into  physical  space  via  the  right  eigenvectors.  The  resulting 
system  of  ODEs  is  advanced  in  time  using  the  third  order  Total  Variation  Diminishing  Runge-Kutta  scheme 
(RK-TVD).  See  [1]  for  further  details  of  the  algorithm. 


4.2.1  Shock-density  wave  interaction 


Consider  the  one  dimensional  Mach  3  shock-entropy  wave  interaction  [5], 
conditions: 


f  (  3.857143,  2.629369,  10.33333  ) 

\  (  1 -I- 0.2sin(fca;),  0,  1  ) 


specified  by  the  following  initial 


— 5  <  a;  <  — 4 
— 4  <  a;  <  5 


(32) 


where  x  S  [—5,  5]  and  k  =  5.  The  solution  of  this  problem  consists  of  a  number  of  shocklets  and  fine  scales 
structures  which  are  located  behind  a  right-going  main  shock. 

Figure  7(a)- (b)  provides  a  comparison  for  all  schemes  at  t  =  1.8  with  an  increasing  number  of  points. 
We  shall  refer  to  the  solution  computed  by  the  WENO-M  scheme  with  N  =  2000  points  as  the  ’’exact” 
solution.  At  a  low  resolution,  N  =  200,  as  shown  in  Eig.  7(a),  WENO-M  and  WENO-Z  capture  much  more 
features  of  the  solution  than  WENO-JS,  particularly  at  the  high  frequency  waves  behind  the  shock,  with 
WENO-Z  achieving  deeper  valleys  and  higher  pikes  than  WENO-M.  Increasing  the  resolution  to  A^  =  300, 
as  shown  in  Fig.  7(b),  we  see  that  convergence  of  WENO-M  and  WENO-Z  is  faster  than  WENO-JS. 

At  Eig.  8,  the  wave  number  k  is  increased  to  10,  yielding  rougher  numerical  approximations  at  the  sine 
wave  density  held  perturbation.  The  same  observation  can  be  made  similar  to  one  made  in  the  previous 
k  =  5  case,  where  the  improvement  of  the  more  accurate  schemes  WENO-M  and  WENO-Z  over  the  classical 
WENO-JS  is  more  distinct. 
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(a)  N  =  200 


(b)  N  =  300 


Figure  7:  Solution  of  the  Mach  3  shock  density  wave  interaction  with  fc  =  5  as  computed  by  WENO-JS, 
WENO-M  and  WENO-Z  schemes  with  (a)  N  =  200,  (b)  N  =  300  points.  The  ’’exact”  solution  is  computed 
by  the  WENO-M  scheme  with  N  =  2000  points. 


Figure  8:  Solution  of  the  Mach  3  shock  density  wave  interaction  with  fc  =  10  computed  by  WENO-JS, 
WENO-M  and  WENO-Z  with  N  =  510  points.  For  clarity,  only  symbols  at  every  fourth  data  point  are 
plotted.  The  ’’exact”  solution  is  computed  by  the  WENO-M  scheme  with  N  =  2000  points. 
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4.2.2  Interacting  Blastwaves 

The  one  dimensional  Blast  waves  interaction  problem  by  Woodward  and  Collela  [7]  has  the  following  initial 


configuration,  with  reflective  boundary  conditions: 

r  ( 

1,  0, 

1000  ; 

)  0  <  X  <  0.1 

{p,U,P)={  ( 

1,  0, 

0.01  ; 

)  0.1<x<0.9  . 

(33) 

1  ( 

1,  0, 

100  ; 

)  0.9<x<1.0 

The  initial  pressure  gradients  generate  two  density  shock  waves  that  collide  and  interact  later  in  time, 
forming  a  profile  as  shown  in  Fig.  9  at  t  =  0.038.  All  three  schemes  converge,  as  the  number  of  points 
increase,  to  the  reference  solution  computed  by  the  WENO-M  with  N  =  4000  points.  As  before,  WENO-M 
and  WENO-Z  show  an  improved  convergence  with  respect  to  WENO-JS,  due  to  their  smaller  dissipation. 
Figure  10  presents  a  separate,  and  more  detailed,  comparison  between  WENO-M  and  WENO-Z,  at  two 
different  portions  of  the  domain  and  different  numbers  of  points.  Careful  examination  of  Fig.  10(a)  shows 
that  WENO-Z  obtains  a  sharper  resolution  of  the  shock  near  x  =  0.65  and  a  higher  pike  near  x  =  0.78.  In 
Fig.  10(b),  the  high  gradient  structure  at  a;  =  0.86  is  clearly  better  resolved  with  WENO-Z,  as  well  as  the 
contact  discontinuity  near  x  =  0.8,  with  N  =  800  points. 


Figure  9:  Solution  of  the  interactive  blastwaves  problem  computed  by  WENO-JS,  WENO-M  and  WENO-Z 
with  N  =  400  points.  For  clarity,  only  symbols  at  every  other  data  point  are  plotted.  The  ’’exact”  solution 
is  computed  by  the  WENO-M  scheme  with  N  =  4000  points. 


5  Conclusions 

We  have  devised  an  improved  version  of  the  fifth  order  WENO  scheme  that  makes  use  of  higher  order 
information  on  the  regularity  of  the  solution  already  contained  in  the  original  framework  of  the  classical 
WENO  scheme.  A  set  of  new  smoothness  indicators  was  devised,  resulting  in  a  less  dissipative  approximation 
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(a)  X  e  [0.625,0.800] 


(b)  a;  G  [0.790, 0.875] 


Figure  10:  Zoom  in  regions  of  the  solution  of  the  interactive  blastwaves  problem  computed  by  WENO-M 
and  WENO-Z  with  N  =  800  points.  For  clarity,  only  symbols  at  every  other  data  point  are  plotted.  The 
’’exact”  solution  is  computed  by  the  WENO-M  scheme  with  N  =  4000  points. 


near  discontinuities.  To  demonstrate  this,  we  have  compared  the  new  WENO  scheme  with  the  classical 
WENO  scheme  of  Jiang  and  Shu  and  the  mapped  WENO  scheme  of  Hendrick  et  al.  in  the  one  dimensional 
solution  of  the  Euler  equations,  in  particular,  the  Mach  3  shock  density  wave  interaction  and  the  interactive 
blastwaves  problems.  These  results  showed  that  the  new  WENO  scheme  substantially  improves  on  the 
classical  WENO  scheme,  generating  numerical  solutions  that  are  slightly  better  than  the  mapped  WENO. 
The  WENO  scheme  proposed  in  this  article  is  a  simple  and  natural  improvement  on  the  classical  smoothness 
indicators  with  no  additional  computational  cost.  Nevertheless,  it  inherits  the  same  theoretical  deficiency 
of  the  classical  WENO  at  critical  points,  no  relevant  influence  at  problems  with  shocks  have  been  found. 
Research  is  currently  underway  to  extend  the  idea  to  higher  order  WENO  reconstruction  schemes  and  will 
be  reported  on  an  upcoming  paper. 
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